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Abstract 

LISA (Laser Interferometer Space Antenna) is a proposed space mission, which will use coherent 
laser beams exchanged between three remote spacecraft to detect and study low-frequency cosmic 
gravitational radiation. In the low-part of its frequency band, the LISA strain sensitivity will be 
dominated by the incoherent superposition of hundreds of millions of gravitational wave signals 
radiated by inspiraling white-dwarf binaries present in our own galaxy. In order to estimate the 
magnitude of the LISA response to this background, we have simulated a synthesized population 
that recently appeared in the literature. Our approach relies on entirely analytic expressions of the 
LISA Time-Delay Interferometric responses to the gravitational radiation emitted by such systems, 
which allows us to implement a computationally efficient and accurate simulation of the background 
in the LISA data. We find the amplitude of the galactic white-dwarf binary background in the 
LISA data to be modulated in time, reaching a minimum equal to about twice that of the LISA 
noise for a period of about two months around the time when the Sun-LISA direction is roughly 
oriented towards the Autumn equinox. This suggests that, during this time period, LISA could 
search for other gravitational wave signals incoming from directions that are away from the galactic 
plane. Since the galactic white-dwarfs background will be observed by LISA not as a stationary but 
rather as a cyclostationary random process with a period of one year, we summarize the theory of 
cyclostationary random processes, present the corresponding generalized spectral method needed to 
characterize such process, and make a comparison between our analytic results and those obtained 
by applying our method to the simulated data. We find that, by measuring the generalized spectral 
components of the white-dwarf background, LISA will be able to infer properties of the distribution 
of the white-dwarfs binary systems present in our Galaxy. 
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I. INTRODUCTION 



The Laser Interferometric Space Antenna (LISA) is a space mission jointly proposed to the 
National Aeronautics and Space Administration (NASA) and the European Space Agency 
(ESA). Its aim is to detect and study gravitational waves (GW) in the millihertz frequency 
band. It will use coherent laser beams exchanged between three identical spacecraft forming 
a giant (almost) equilateral triangle of side 5 x 10 6 kilometers. By monitoring the relative 
phase changes of the light beams exchanged between the spacecraft, it will extract the 
information about the gravitational waves it will observe at unprecedented sensitivities l|. 

The astrophysical sources that LISA is expected to observe within its operational fre- 
quency band (10~ 4 — 1 Hz) include extra-galactic super-massive black-hole coalescing bi- 
naries, stochastic gravitational wave background from the early universe, and galactic and 
extra-galactic coalescing binary systems containing white dwarfs and neutron stars. 

Recent surveys have uniquely identified twenty binary systems emitting gravitational 
radiation within the LISA band, while population studies have concluded that the large 
number of binaries present in our own galaxy should produce a stochastic background that 
will lie significantly above the LISA instrumental noise in the low-part of its frequency 
band. It has been shown in the literature (see 3] for a recent study and 0, 0] for earlier 
investigations) that these sources will be dominated by detached white-dwarf — white-dwarf 
(WD- WD) binaries, with 1.1 x 10 8 of such systems in our Galaxy. The detached WD- WD 
binaries evolve by gravitational-radiation reaction and the number of such sources rapidly 
decreases with increasing orbital frequency. Although it is expected that, above a certain 
frequency cut-off (1 — 2 mHz), we will be able to resolve individual signals and remove 
them from the LISA data, it is still not clear how to further improve the LISA sensitivity 
to other gravitational wave signals in the region of the frequency band below the WD- WD 
background frequency cut-off. Although two promising data analysis procedures have been 
proposed 0, ||| for attempting to subtract the galactic background, considerable work still 
needs to be done to verify their effectiveness. In this context, simulating the LISA response 
to the WD- WD background will be particularly useful for verifying present and future data 
analysis "cleaning" algorithms. A realistic simulation will also quantify the effects of the 
LISA motion around the Sun on the overall amplitude and phase of the GW signal generated 
by the background in the LISA data. The directional properties of the LISA response and 
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its time dependence introduced by the motion of LISA around the Sun, together with the 
non-isotropic and non-homogeneous distribution of the WD- WD binary systems within the 
galactic disk as seen by LISA, imply that the magnitude of the background observed by 
LISA will not be a stationary random process. As a consequence of the one- year periodicity 
of the LISA motion around the Sun, there exist relatively long (« 2 months) stretches of 
data during which the magnitude of the LISA response to the background will reach an 
absolute minimum |7]. Our simulation shows this minimum to be less than a factor of two 
larger than the level identified by the LISA secondary noises, suggesting the possibility of 
performing searches for gravitational radiation from other sources located in regions of the 
sky that are away from the galactic plane. The LISA sensitivity to such signal in fact will 
be less limited by the WD- WD background during these periods of observation. 

This paper is organized as follows. In Section [H] we provide the analytic expression of one 
of the LISA Time-Delay Interferometric (TDI) responses to a signal radiated by a binary 
system. Although all the TDI responses to binary signals were first given in their closed 
analytic form in !$], in what follows we will focus our attention only on the unequal-arm 
Michelson combination, X. In Section ITTTI we give a summary of how the WD- WD binary 
population was obtained, and a description of our numerical simulation of the X response 
to it. In Section IIVI we describe the numerical implementation of our simulation of the 
LISA X response to the WD- WD background, and summarize our results. In particular, 
in agreement with the results by Seto 7], we find the amplitude of the galactic WD- WD 
background in the LISA X-combination to be modulated in time, reaching a minimum when 
the Sun-LISA direction is roughly oriented towards the Autumn equinox. Furthermore, we 
show the amplitude of the background at its minimum to be a factor less than two larger than 
the level identified by the LISA noise for a time period of about two months, suggesting that 
LISA could search (during this time period) for other gravitational wave signals incoming 
from regions of the sky that are away from the galactic plane. 

The time-dependence and periodicity of the magnitude of the WD- WD galactic back- 
ground in the LISA data implies that it is not a stationary but rather a cyclostationary 
random process of period one year. After providing a brief summary of the theory of cy- 
clostationary random processes relevant to the LISA detection of the WD- WD galactic 
background, we apply it to three years worth of simulated LISA X data. We find that, 
by measuring the generalized spectral components of such cyclostationary random process, 
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FIG. 1: Schematic LISA configuration. Each spacecraft is equidistant from point o, with unit 
vectors pi indicating directions to the three spacecraft. Unit vectors hi point between spacecraft 
pairs with the indicated orientation. 

LISA will be able to infer key-properties of the distribution of the WD- WD binary systems 
present in our own Galaxy. 



II. THE LISA RESPONSE TO SIGNALS FROM BINARY SYSTEMS 

The overall LISA geometry is shown in Figure There are six beams exchanged be- 
tween the LISA spacecraft, with the six Doppler measurements = 1,2,3) recorded 
when each received beam is mixed with the laser light of the receiving optical bench. The 
frequency fluctuations from the six lasers, which enter in each of the six Doppler measure- 
ments, need to be suppressed to a level smaller than that identified by the secondary (proof 
mass and optical path) noises [t| in order to detect and study gravitational radiation at the 
predicted amplitudes. 

Since the LISA triangular array has systematic motions, the two one-way light times 
between any spacecraft pair are not the same 10]. Delay times for light travel between the 
spacecraft must be accounted for depending on the sense of light propagation along each 
link when combining these data as a consequence of the rotation of the array. Following 
the arms are labeled with single numbers given by the opposite spacecraft; e.g., arm 2 (or 
2') is opposite spacecraft 2, where primed delays are used to distinguish light-times taken 
in the counter-clockwise sense and unprimed delays for the clockwise light times (see Figure 
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FIG. 2: Schematic diagram of LISA configurations involving six laser beams. Optical path delays 
taken in the counter-clockwise sense are denoted with a prime, while unprimed delays are in the 
clockwise sense. 

(J2J)). Also the following labeling convention of the Doppler data will be used. Explicitly: 
1/23 is the one-way Doppler shift measured at spacecraft 3, coming from spacecraft 2, along 
arm 1. Similarly, y^ is the Doppler shift measured on arrival at spacecraft 2 along arm 1' 
of a signal transmitted from spacecraft 3. Due to the relative motion between spacecraft, 
Li 7^ L\ in general. As in 0, 12 1, we denote six further data streams, = 1,2,3), as 



the intra-spacecraft metrology data used to monitor the motion of the two optical benches 
and the relative phase fluctuations of the two lasers on each of the three spacecraft. The 
frequency fluctuations introduced by the lasers, by the optical benches, by the proof masses, 
by the fiber optics, and by the measurements themselves at the photo-detectors (i.e. the 
shot-noise fluctuations) enter the Doppler observables yij, with specific time signatures; 



see Refs. (3, [l2j for a detailed discussion. The contribution y^ w due to GW signals was 
derived in Ref. in the case of a stationary array, and further extended to the realistic 
configuration js| of the LISA array orbiting around the Sun. 

Let us consider for instance the "second generation" unequal-arm Michelson TDI observ- 
ables, (Xi,X 2 ,X 3 ). Their expressions, in terms of the Doppler measurements y^, Zy, are as 
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follows [1 



Xi = [{y 3 i + yn-2) + (2/21 + yi2;3') ;2 '2 + (^ 21 + 2/l2;3');33'2'2 ~^ (^ 31 2/ 13 ; 2 )-33'33'2'2] 
-[(l/21 + J/l2;3') + (2/31 + £/l3;2). 3 3, + (2/31 + 2/l3;2) 2 '233' + (2/21 + 2/l2;3');2'22'233'] 
+ 2 [( Z 21 - z 3l) - (^21 - ^3l) ; 33' - (^21 ~ ^3l) ;2 '2 + ( Z 21 ~ 2 3l) ;3 3' 3 3'2'2 
+ (Z 2 1 - ^3l) ;2 '22'233' ~~ ( Z 21 ~ Z 3l) ;2 '233'33'2'2] ' W 

with X 2 , X 3 following from Eq. (Q) by permutations of the spacecraft indices. The semicolon 
notation shown in equation Q emphasizes that the operation of sequentially applying two 
or more delays to a given measurement is non-commutative as consequence of the time 
dependence of the light-times Li and L\ (i = 1, 2, 3), and a specific order has to be adopted 
to adequately suppress the laser noises (ll|, |l5( . Specifically: y^u = Vij(t — Li{t) —L k (t — 
Li)) 7^ yij-ik (units in which the speed of light c = 1). 

The expressions of the gravitational wave signal and the secondary noise sources entering 
into Xi will in general be different from those entering into X, the corresponding "first 
generation" unequal-arm Michelson observable derived under the assumption of a stationary 
LISA array Q, However, the magnitude of the corrections introduced by the motion of 
the array are proportional to the product between the time derivative of the GW amplitude 
and the difference between the actual light travel times and those valid for a stationary array. 
At 1 Hz, for instance, the larger correction to the signal (due to the difference between the co- 
rotating and counter-rotating light travel times) is two orders of magnitude smaller than the 
main signal. Since the amplitude of this correction scales linearly with the Fourier frequency, 
we can completely disregard this effect (and the weaker effect due to the time dependence 



of the light travel times) over the entire LISA band n| . Furthermore, since along the LISA 
orbit the three armlengths will differ at most by ~ 1 %-2%, the degradation in signal-to- noise 
ratio introduced by adopting signal templates that neglect the inequality of the armlengths 
will be of only a few percent. For these reasons, in what follows we will focus on the 
expressions of the GW responses of various second-generation Time-Delay Interferometry 
(TDI) observables by disregarding the differences in the delay times experienced by light 
propagating clockwise and counterclockwise, and bv assuming the three LISA armlengths to 
be constant and equal to L = 5 x 10 6 km ~ 16.67 s These approximations, together with 
the treatment of the moving-LISA GW response discussed in |8( are essentially_^equivalent 
to the rigid adiabatic approximation of Ref. 



)onse discussed in |8( are essentially e 
[l^ , and to the formalism of Ref. p| • 



These considerations imply that the second generation TDI expressions for the gravita- 
tional wave signal and the secondary noises can be expressed in terms of the corresponding 
first generation TDIs. For instance, the gravitational wave signal entering into the sec- 
ond generation unequal-arm Michelson combination, Xp w , can be written in terms of the 
gravitational wave response of the corresponding first generation unequal-arm Michelson 
combination, X GW (t), in the following manner |l7j | 

X GW (t) = X GW {t) - X GW {t - 4L) (2) 

Equation (J2J) implies that any data analysis procedure and algorithm that will be imple- 
mented for the second generation TDI combinations can actually be derived by considering 
the corresponding first generation TDI expressions. For this reason, from now on we will 
focus our attention on the gravitational wave responses of the first generation combinations. 

The gravitational wave response X GW (t) of the unequal-arm Michelson TDI combination 
to a signal from a binary system has been derived in and it can be written in the following 
form 

X GW (t) = ft[A(x,t) e-**<*>] , (3) 

where x = u s L (u s being the angular frequency of the GW signal in the source reference 
frame), and the expressions for the complex amplitude A(x,t) and the real phase (f)(t) are 

A(x,t) = 2x sin(x) { [smc[(l + c 2 (t))|] e ix ^ +d2(t)) + sinc{{\ - c 2 (t))|] e ix ^ +d ^ B 2 (t) 
sinc[(l - c 3 (t))|] e^i+^M) + sinc[(l + c 3 (t))|] e^i+^J B 3 (t)} , (4) 



<p(t) = oj s t + oj s R cos (3 cos(u; s t + rjo — A) . (5) 

In equation (0) R is the distance of the guiding center of the LISA array, o, from the Solar 
System Barycenter, (/3, A) are the ecliptic latitude and longitude respectively of the source 
location in the sky, Q = 27r/year, and 7] defines the position of the LISA guiding center in 
the ecliptic plane at time t — 0. Note that the functions Cfc(t), dk(t), and Bk(t) (k = 2, 3) do 
not depend on x. The analytic expressions for Ck{t), and dk(t) are the same as those given 
in equations (46,47) of reference |8j, while the functions Bk(t) (k = 2,3) are equal to 

B k (t) = (a^ + % a (3) ) u k (t) + (a^ + i a^) v k (t) . (6) 



S 



The coefficients (a^\ a^ 2 \ a^ 3 \ a' 4 ') depend only on the two independent amplitudes of the 
gravitational wave signal, (h + , h x ), the polarization angle, ip, and an arbitrary phase, <fio, 
that the signal has at time t = 0. Their analytic expressions are given in equations (41-44) 
of reference U, while the functions Uk(t), and (k = 2, 3) are given in equations (27,28) 
in the same reference. 

Since most of the gravitational wave energy radiated by the galactic WD- WD binaries will 
be present in the lower part of the LISA sensitivity frequency band, say between 10~ 4 — 10~ 3 
Hz, it is useful to provide an expression for the Taylor expansion of the X response in the 
long- wavelength limit (LWL), i.e. when the wavelength of the gravitational wave signal 
is much larger than the LISA armlength (x « 1). As it will be shown in the following 
sections, the LWL expression will allow us to analytically describe the general features of 
the white dwarfs background in the X-combination, and derive computationally efficient 
algorithms for numerically simulating the WD- WD background in the LISA data. 

The nth-order truncation, X9W(t), of the Taylor expansion of X GW (t) in power series of 
x can be written in the following form 

n 

XgW(t) = Re A (k) (t) x k+2 e-^W , (7) 

k=0 

where the first three functions of time A^ k \t), k < 2 are equal to 

A® = 4 [B 2 - B 3 \ , 

A® = At [(d 2 + 2)B 2 - (d 3 + 2)B 3 ] , 

A^ = [2d, 2 + 8d 3 + ^ + \c 3 2 ] B 3 - [2d 2 2 + 8d 2 + ^ + \c 2 2 } B 2 . (8) 

Note that the form we adopted for X GW (t) (equation^ makes the derivation of the functions 
A^ k \t) particularly easy since the dependence on x in A(x,t) is now limited only to the 
coefficients in front of the two functions B 2 (t) and B 3 {t) (see equation (jlj)). 

Although it is generally believed that the lowest order long-wavelength expansion of the X 
combination, X^, is sufficiently accurate in representing a gravitational wave signal in the 
low-part of the LISA frequency band, there has not been in the literature any quantitative 
analysis of the error introduced by relying on such a zero-order approximation. Since any 
TDI combination will contain a linear superposition of tens of millions of signals, it is crucial 
to estimate such an error as a function of the order of the approximation, n. In order to 
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FIG. 3: Plots of the percentage root-mean-squared errors, M{X GW ,X™), associated with the 
long- wavelength expansion index n, as functions of the gravitational wave frequency, f s . The source 
location has been assumed to be in the center of our galaxy. 

determine how many terms we need to use for a given signal angular frequency, u 8 , we will 
rely on the following " matching function" 



M(X 



GW Y GW\ 



^[X^(t)-X^(t)] 2 dt 



(9) 



fi[XW(t)) 2 dt 

Equation estimates the percent root-mean-squared error implied by using the n th order 
LWL approximation. In Figure (jHJ) we plot M as a function of the signal frequency, f s 
(= lj s /2tt), for n = 0, 1, 2. At 5 x 10 -4 Hz, for instance, the zero-order LWL approximation 
(n = 0) of the X combination shows an r.m.s. deviation from the exact response equal to 
about 10 percent. As expected, this inaccuracy increases for signals of higher frequencies, 
becoming equal to 40 percent at 2 x 10~ 3 Hz. With n — 1 the accuracy improves showing 
that the Xffr response deviates from the exact one with an r.m.s. error smaller than 10 
percent in the frequency band (10 -4 — 2 x 10 -3 ) Hz. In our simulation we have actually 
implemented the n = 2 LWL expansion because it was possible and easy to do. 
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III. WHITE DWARF BINARY POPULATION DISTRIBUTION 



The gravitational wave signal radiated by a WD- WD binary system depends on eight 
parameters, (<f) , t, ip, D, (3, A, M. c , u s ), which are the constant phase of the signal (4> ) at 
the starting time of the observation, the inclination angle (l) of the angular momentum 
of the binary system relative to the line of sight, the polarization angle (ip) describing the 
orientation of the wave polarization axes, the distance (D) to the binary, the angles (X,/3) 
describing the location of the source in the sky relative to the ecliptic plane, the chirp mass 
(Ai c ), and the angular frequency (uj s ) in the source reference frame respectively. Since it can 
safely be assumed that the chirp mass Ai c and the angular frequency u s are independent of 
the source location |2j and of the remaining angular parameters <j) , i, ip, and because there 
are no physical arguments for preferred values of the constant phase <p and the orientation 
of the binary given by the angles l and ip, it follows that the joint probability distribution, 
P(4> , i>i D, (3, A, M. c , u s ), can be rewritten in the following form 

P{(j) ,L,^,D,(5,\,M c ,u s ) = P 1 { ( f) )P 2 {L)P 3 {^)P i {D,l3,X)P 5 {M c ,u s ) . (10) 

In the implementation of our simulation we have assumed the angles O and ip to be uniformly 
distributed in the interval [0, 2ir), and cost uniformly distributed in the interval [—1, 1]. We 
further assumed the binary systems to be randomly distributed in the Galactic disc according 
to the following axially symmetric distribution V/±{R,z) (see ^ Eq. (5)) 

V( P A- e ~ R/H sech\z/z ) 

where (R, z) are cylindrical coordinates with origin at the galactic center, H = 2.5 kpc, and 
z Q = 200 pc, and it is proportional to P A (D,\,(3) through the Jacobian of the coordinate 
transformation. Note that the position of the Sun in this coordinate system is given by 
Rq = 8.5 kpc and z & = —30 pc. We then generate the positions of the sources from the 
distribution given by Eq. (fTT|) and map them to their corresponding ecliptic coordinates 

(AAA). 

The physical properties of the WD- WD population (M. c = (mim2) 3 ' 5 / (mi + m2) 1//5 , with 
mi, m 2 being the masses of the two stars, and u s = 2nf s = 47r/orbitalperiod) are taken from 
the binary population synthesis simulation discussed in Il8l . For details on this simulation 
we refer the reader to , and for earlier work to P, [J 0, lioi 2ol I21I ] . The basic ingredient for 
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these simulations is an approximate binary evolution code. A representation of the complete 
Galactic population of binaries is produced by evolving a large (typically 10 6 ) number of 
binaries from their formation to the current time, where the distributions of the masses 
and separations of the initial binaries are estimated from the observed properties of local 
binaries. This initial-to-final parameter mapping is then convolved with an estimate of the 
binary formation rate in the history of the Galaxy to obtain the total Galactic population 
of binaries at the present time. From these the binaries of interest can then be selected. In 
principle this technique is very powerful, although the results can be limited by the limited 
knowledge we have on many aspects of binary evolution. For WD- WD binaries, the situation 
is better than for many other populations, since the observed population of WD- WD binaries 
allows us to gauge the models (e.g. P). 

We also include the population of semi-detached WD- WD binaries (usually referred to 
as AM CVn systems) that are discussed in detail in Q|. In these binaries one white dwarf 
transfers its outer layers onto a companion white dwarf. Due to the redistribution of mass 
in the system, the orbital period of these binaries increases in time, even though the angular 
momentum of the binary orbit still decreases due to gravitational wave losses. The formation 
of these systems is very uncertain, mainly due to questions concerning the stability of the 
mass transfer (e.g. 2^|) 

From the models of the Galactic population of the detached WD- WD binaries and AM 
CVn systems two dimensional histograms were created, giving the expected number of both 
WD- WD binaries and AM CVn systems currently present in the Galaxy as function of the 
log of the GW radiation frequency, / s (= ^s/27r) and chirp mass, M. c . In the case of the 
detached WD- WD binaries, the (log/ s , M. c )) space was defined over the set M. c G (0, 1.5], 
log/ s G [—6,-1], and contained 30 x 50 grid points, while in the case of the AM CVn 
systems the region is intrinsically smaller, M. c G (0, 1.2], log/ s G [—4, —1.5], containing only 
24 x 25 grid points. 

Figure (jlj) shows the distribution of the number of detached WD-WD binaries as a func- 
tion of the chirp mass and signal frequency in the form of a contour plot. This distribution 
reaches its maximum within the LISA frequency band when the chirp mass is equal to 
~ 0.25 A^0, and it monotonically decreases as a function of the signal frequency. The dis- 
tribution of the number of AM CVn systems has instead a rather different shape, as shown 
by the contour plot given in Figure (0). The region of the (Ai c , log/ s ) space over which the 
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FIG. 4: The distribution of detached white-dwarf — white-dwarf binaries in our galaxy as a 
function of the gravitational wave frequency, f s , and chirp mass, .M re- 
distribution is non-zero is equal to (0,0.07) x (—3.4, —2.2), and it reaches its maximum at 
the point (0.03,-3.35). 

IV. SIMULATION OF THE BACKGROUND SIGNAL IN THE LISA DATA 

In order to simulate the LISA X response to the population of WD-WD binaries derived 
in Section IHII one needs to coherently add the LISA response to each individual signal. 
Although this could naturally be done in the time domain, the actual CPU time required 
to successfully perform such a simulation would be unacceptably long. The generation in 
the time domain of one year of X GW (t) response to a single signal sampled at a rate of 16 
seconds would require about 1 second with an optimized C++ code running on a Pentium 
IV 3.2 GHz processor. Since the number of signals from the background is of the order 10 8 , 
it is clear that a different algorithm is needed for simulating the background in the LISA 
data within a reasonable amount of time. We were able to derive and implement numerically 
an analytic formula of the Fourier transform of each binary signal, which has allowed us to 
reduce the computational time by almost a factor 100. Furthermore, we have run our code 
on the Jet Propulsion Laboratory (JPL) supercomputer system, which includes 64 Intel 
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FIG. 5: The distribution of AM CVn binary systems in our galaxy as a function of the gravitational 
wave frequency, f s , and chirp mass, M. c . 

Itanium2 processors each with a clock speed of 900 MHz. 



A. The Fourier transform of a binary signal 

The expression of the Fourier transform of the TDI response X GW to a single binary signal 
(Eqs. (0 IH EI), cannot be written (to our knowledge) in closed analytic form. However, 
by using the LWL expansion of the X GW -response, it is possible to obtain a closed-form 
expression of its Fourier transform. Since the WD- WD binary background has a natural 
frequency cut-off that is between 1 and 2 millihertz, the LWL expansion of the X GW response 
(EqEJ), truncated at n = 2, can be used for accurately representing the gravitational wave 
response of each binary signal, as discussed in Section |Hj 

In order to derive the Fourier transform of X%^(t), we use the following expansion of 
the function e~ 1 ^ in terms of the Bessel functions of the first kind, J q , j3| 

oo 

e -i 4>(t) = J2 J q {u s Rcosf3) e- i[uJst + fl (o* + w - A + f)] _ ( 12 ) 

<?=— oo 

Since the Bessel functions | J q (u s R cos (3)\ are much smaller than unity when \q\ >> 
\u s Rcos P\, the expansion given by equation (fT2*|) can be truncated at a finite index Q, 
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providing an accurate numerical estimation of the function e~ l These considerations 
allow us to write the following expression of the Fourier transform of X?^ (t) 

( n Q 

Xffi{u) = T i^Yl Yl A(k) ^ x " +2 J q (usRcoa/3) + i + V o - x + f) 

\ k=0 q=-Q 
Q n 

= 7T J q (oo S R cos (3) xk+2 {^{^[^(t)]) * + u s + q9)e iq ^ 0+x -*W 

q=-Q k=0 

+ 5(-w + lo s + q9)e iq ^~ x+v/2) ] + iT (S[A (fc) (t)]) * [5(u + oo s + qny^+^M 
- 8(-oj + oo s + gfi)e^°- A+7r/2) ] } , (13) 

where JF is the Fourier transform operator, the symbol * between two expressions means 
their convolution, uj is the Fourier angular frequency, and (t) are defined in Eq. (JJJ) and 
given in Eq. © with k — 0, 1, 2. 

hie rronQvol tavvyihIo bw +- li /z* Hr\nnor 4" v o m cmr rvi r\T f ho V 



As an example application of this general formula for the Fourier transform of the Xpy*(t) 



response, let us apply it to the lowest order LWL expansion (n = 0) 

Q 

X^{uj) = 4:71 5Z J q( UJ sRcosP)x 2 {[ai(u 2 (uj) -u 3 (u)) + a 2 (v 2 {u) - v 3 (u))] 

q=-Q 

* [5(u + u s + q{l)e iq{ -' no+x ~ 7T/2) + + uo s + qQ)e iq{v °- x+7r/2) ] 
+i [a 3 (u 2 {u)) - u 3 (u)) + ai(v 2 (uj) - Vs(u))] 

* [5{uj + u s + qQ)e iq{ ~ V0+x - 7T/2) - 5{-uj + uj s + qQ)e iq{r '°- x+w/2) } } . (14) 

Since the Fourier transforms of u and v are both linear combinations of nine Dirac delta 
unctions centered on the frequencies ±/ Q , I = 0, 1, 2, 3, 4 (see equations (27-30) in reference 
^1 for the expressions of u and v), it follows that X^(u) is also a linear combination of Dirac 
delta functions. In particular, in the limit of negligible Doppler modulation, the resulting 
expression (fT4^) reduces as expected to that of a purely amplitude modulated sinusoidal 
signal with central frequency equal to u s and upper and lower band-limits given by u s + 4fl 



and uj s — 4f2 respectively 



24 1 . 



The actual expression of the Fourier transform we implemented in our simulation of the 
WD- WD background used Eq. (|13|) with n = 2, and maximum value of the index of the 
Bessel expansion, Q, equal to \u s R cos/3| +20 in order to make negligible the error associated 
with the truncation of the expansion itself. 

One extra mathematical detail that we need to include is that the Fourier transform 
of X9jW(t) is performed over a finite integration time, T, while the expression in Eq. (|13|) 
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corresponds to an infinite-time Fourier transform. In order to account for this discrepancy we 
convolved the analytic Fourier transform of the signal given in equation (|T3"j) with the Fourier 
transform of a window function with an integration time T. To avoid leakage introduced by 
using a simple rectangular window, we have used instead the Nuttall's modified Blackman- 
Harris window 25]. Although this window is characterized by having the main lobe of its 
Fourier transform slightly wider than that of the rectangular window, the maximum of its 
side-lobes are about four orders of magnitude lower than those of the rectangular window, 
reducing leakage significantly. The expression of its Fourier transform, W(u), is equal to 

W(w) = n [Sinc(ujT) +iCosc(ujT)} (15) 

— ni [Sinc(ujT + 2ix) + Sinc(uT — 2n) + i [CosciujT + 27r) + Cosc(u>T — 2ir))] 

— n 2 [Sinc{uoT + 4vr) + Sinc{uT - 4vr) + % (Cosc{uoT + 4tt) + Cosc{uT - 4tt))] 

— n 3 [SinciuT + 6tt) + SinciuT — 6n) + i (Cosc{ujT + 6n) + Cosc(ujT — 6%))} , 

where the functions Sinc(.) and Cosc(.) are defined as follows 

S inc (.) s , Cosc(.) S C ° s( - } ~ 1 , (16) 

and the coefficients n r , r = 0, 1, 2, 3 have the following numerical values 

n = 0.3635819 , n x = 0.24458875 , n 2 = 0.06829975 , n 3 = 0.00532055 . (17) 



B. Generation of the signal parameters 

We used the distributions given in Section ITTT1 to randomly generate the parameters <p , l, 
ip, D, (3, A, while the values of the chirp mass, M. Cl and the logarithm of the frequency of the 
signal, log(/ s ), were obtained by further processing the numeric distribution function (given 
in Section ITTT|) of the number of sources. To derive the distribution function for the variables 
{M.c-1 logCfs)) within each grid-rectangle of our numerical distribution we proceeded in the 
following way Let us consider the number of sources N(xi,X2) as a function of two 

coordinates (xi,x 2 ) of a point in the {M. c , log(/ s )) plane within a specified grid-rectangle 
of the numerical distribution. This function can be approximated there by the following 
quadratic polynomial 

N(xi, x 2 ) = n o(l - xx)(l - x 2 ) + n u xix 2 + n 01 (l - xi)x 2 + n w xx{l - x 2 ) , (18) 
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where noo, nn, n i, nio are equal to the number of signals at the "corners" of the considered 
grid-rectangle and are obtained by interpolation; (xi,x 2 ) are therefore two real numbers 
defined in the range [0, 1]. 

If we integrate out the ^-dependence in N(x\, x 2 ) we obtain 

Ar f x f 1 ATt SJ (-^oo + riu - riQi + n 1Q )xi + n 00 + n i 

Nxitxi) = J N(xi,x 2 )dx 2 = , (19) 

which defines the total number of sources within that grid-rectangle having chirp mass equal 
to x\. In order to derive the probability distribution function of x\ within that grid-rectangle 
we can define the following mapping between a uniformly distributed random variable, say 
z\, and the random variable X\ 

z _ Ip 1 N X2 (x , 1 )dx , 1 _ (-n 00 + nn - n 01 + n 1Q )x\ + 2(n 00 + ngijxi 
1 Jq 1 N X2 (x' l )dx' l n 00 + nn + n 01 + n 10 

By solving the above non-linear equation for every uniformly sampled Z\ we obtain 

n o + npi - A/ng + 2n 00 n 01 + n 2 Ql + (— ra§ - 2n 00 n 01 + n 2 u + 2n u n w - n 2 Q1 + n 2 10 )z 1 

n 00 - nn + n i - ni 

(21) 

where the branch "-" has been chosen such that X\ remains in the range [0, 1]. If noo — n ii + 
^oi — ^io = then equation (|2~T|) is no longer valid, and we have instead x\ — Z\. 

A similar procedure can be implemented for calculating x 2 - By integrating N(xi,x' 2 ) 
with respect to x' 2 over the range (0, x 2 ), we can establish the following relationship between 
another uniformly distributed random variable, say z 2 , and x 2 

Ip 2 N(x 1 ,x' 2 )dx' 2 _ [(nop + n lx - n i - n 10 )xi - n 00 + n i]x| + 2[(-n 00 + ni )xi + n 0Q }x 2 
Jg 1 N(x u x' 2 )dx' 2 {-noo + n n - n i + ni )xi + n 00 + n i 

(22) 

After some simple algebra we can finally solve for x 2 in terms of x\ (itself a function of the 
uniformly distributed random variable z\) and z 2 



z 2 



_ noo(xi - 1) - n 10 x 1 + ^/F(x!)z 2 + (n 2 m - 2n oni + n\ Q )x\ + 2(-ng + n oni )a;i + n 2 m 

(nn - ^oi - n 10 + n o)x 1 - n 00 + n i 

(23) 

where now we have chosen the "+" branch so x 2 G [0, 1] range, and the function F(xi) is 
equal to 

F(xx) = (n 2 u +nl 1 -nl -2no 1 n 11 +2noon w -n 2 w )xl+2(-noon w +nl Q -nl^^ . 

(24) 
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Note that the equation for x<i above is no longer valid when (nn — noi — nio +^00)^1 — ^00 + 
riQi = 0, in which case 22 = z-i- Once x\ and X2 are calculated, they can be converted into 
the physical parameters A4 C and uj s according to the following relationships 

M c = M c00 + x x A Mc (25) 
oo s = 2irlO log{fsoo)+X2 Alo ^> (26) 

where (Ai c00 , log(/ s00 )) are the coordinates of the "lower-left-hand-corner" of the considered 
grid-rectangle, and A Mc , A log (j s ) are the lengths of the sides of the grid-rectangle. 



C. Results of the Numerical Simulation 

The expression for the finite-time Fourier transform of each WD-WD signal in the X GW 
response given in Section II V Al allows us to coherently add in the Fourier domain all the 
signals radiated by the WD-WD galactic binary population described in Section 11111 After 
inverse Fourier transforming the synthesized response and removing the window from it, we 
finally obtain the time-domain representation of the background as it will be seen in the 
LISA TDI combination X. This is shown in Figure (|7J), where we plot three years worth 
of simulated X GW (t), and include the LISA noise l|. The one-year periodicity induced 
by the motion of LISA around the Sun is clearly noticeable. One other interesting feature 
shown by Figure (J7j) is that the amplitude response reaches absolute minima when the Sun- 
LISA direction is roughly oriented towards the Autumn equinox, while the absolute maxima 
take place when the Sun-LISA direction is oriented roughly towards the Galactic center [t| . 
This fact can easily be understood by looking at Figure ©• Since the ecliptic plane is not 
parallel to the galactic plane, and because our own solar system is about 8.5 kpc away from 
the galactic center (where most of the of WD-WD binaries are concentrated ) it follows that 
the LISA X GW response does not have a six-months periodicity. 

Note also that, for a time period of about 2 months, the absolute minima reached by 
the amplitude of the LISA response to the WD-WD background is only a factor less than 2 
larger than the level of the instrumental noise. This implies that during these observation 
times LISA should be able to search for other sources of gravitational radiation that are 
not located in the galactic plane. This might turn out to be the easiest way to mitigate 
the detrimental effects of the WD-WD background when searching for other sources of 
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LISA Normal 



to the Galactic Center 




LISA Normal 

FIG. 6: Two snapshots of LISA along its trajectory as recorded six months apart by an observer 
in the ecliptic plane. They correspond to the times when the LISA response to the background 
achieves a local maximum. The magnitudes of these maxima are not equal due to the relative 
disposition of the ecliptic plane with respect to the galactic plane. At time to the angle 9q between 
the normal to the plane of LISA and a vector pointing to the galactic center is equal to 24.47°. Six 
months later, when the LISA response is also at a maximum, the angle 9\ is equal to 35.53° which 
results in a smaller maximum. 

gravitational radiation. We will quantitatively analyze in a follow up work how to take 
advantage of this observation in order to optimally search, during these time periods, for 
sources that are off the galactic plane. 

In Figure JSj) we plot, as functions of the Fourier frequency, /, the windowed Fourier 
powers of both the signal and the noise entering into the TDI X combination. Note that in 
the region of the LISA band below 0.2 millihertz the power of the WD- WD background is 
smaller than that of the instrumental noise. 



V. CYCLOSTATIONARY PROCESSES 

The results of our simulation (Figure (J7J)) indicate that the LISA X GW response to the 
background can be regarded, in a statistical sense, as a periodic function of time. This is 
consequence of the deterministic (and periodic) motion of the LISA array around the Sun. 
Since its autocorrelation function will also be a periodic function of period one year, it fol- 
lows that any LISA response to the WD- WD background should no longer be treated as a 
stationary random process but rather as a periodically correlated random process. These 
kind of processes have been studied for many years, and are usually referred to as cyclo- 



stationary random processes (see 



for a comprehensive overview of the subject and for 



more references). In what follows we will briefly summarize the properties of cyclostationary 
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FIG. 7: Three years of simulated TDI X response to the WD- WD Galactic background signal. 
The time series of the LISA instrumental noise is displayed for comparison. 

processes that are relevant to our problem. 

A continuous stochastic process X(t) having finite second order moments is said to be 
cy do stationary with period T if the following expectation values 



E[X{t)\ = m{t) = m(t + T), 
E[X{t')X(t)\ = C{t',t) = C{? + T,t + T) 



(27) 
(28) 



are periodic functions of period T, for every (£', t) 6 R x R. For simplicity from now on we 
will assume m(t) = 0. 

If X(t) is cyclostationary, then the function B(t,r) = C(t + r,t) for a given r G R is 
periodic with period T, and it can be represented by the following Fourier series 



B(t,r) = J2 B r{ry 2 ^ 



(29) 
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Background and LISH noise in X 
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FIG. 8: The amplitude of the Fourier transform of the WD- WD Galactic background gravitational 
wave signal and of the LISA instrumental noise entering into the TDI combination X. 



where the functions B r (r) are given by 



B r (r) 



T 



B(t,T)e~ l27Tr T dt 



(30) 



The Fourier transforms g r (f) of B r (r) are the so called "cyclic spectra" of the cyclostationary 
process X{t) [27] 

/CO 
B r (r)e- i2 ^ dr . (31) 
-oo 

If a cyclostationary process is real, the following relationships between the cyclic spectra 
hold 



B_ r (r) = B*{t) 

g-r(-f) = <£(/) 



(32) 
(33) 



where the symbol * means complex conjugation. This implies that, for a real cyclostationary 
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process, the cyclic spectra with r > contain all the information needed to characterize the 
process itself. 

The function a 2 {r) = B(0, r) is the variance of the cyclostationary process X(t), and it 
can be written as a Fourier decomposition as a consequence of Eq. (|3U|) 



a 2 



» = Yl R r^i (34) 



where if r = -B r (0) are harmonics of the variance a 2 . From Eq. it follows that H_ T = H*. 

For a discrete, finite, real time series X t , t = 1, . . . , M we can estimate the cyclic spectra 
by generalizing standard methods of spectrum estimation used with stationary processes. 
Assuming again the mean value of the time series X t to be zero, the cyclic autocorrelation 
sequences are defined as 



It has been shown |27( that the cyclic autocorrelations are asymptotically (i.e. for N — > oo) 
unbiased estimators of the functions B r (r). The Fourier transforms of the cyclic auto- 
correlation sequences s[ are estimators of the cyclic spectra g r {f)- These estimators are 
asymptotically unbiased, and are called "inconsistent estimators" of the cyclic spectra, i.e. 
their variances do not tend to zero asymptotically. In the case of Gaussian processes 
consistent estimators can be obtained by first applying a lag window to the cyclic autocor- 
relation and then perform a Fourier transform. This procedure represents a generalization 



3- 



of the well-known technique for estimating the spectra of stationary random processes 

An alternative procedure for identifying consistent estimators of the cyclic spectra is to 
first take the Fourier transform, X(f), of the time series X(t) 

X(f) = J2 X te- i2nf{t - 1} (36) 
t=i 

and then estimate the cyclic periodograms g r (f) 

X{f)X*{f- 2 -^) 
9r(f) = — • (37) 

By finally smoothing the cyclic periodograms, consistent estimators of the spectra g r (f) are 
then obtained. The estimators of the harmonics H r of the variance a 2 of a cyclostationary 
random process can be obtained by first forming a sample variance of the time series X t . 
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The sample variance is obtained by dividing the time series X t into contiguous segments of 
length r such that r is much smaller than the period T of the cyclostationary process, and 
by calculating the variance a 2 over each segment. Estimators of the harmonics are obtained 
either by Fourier analyzing the series o 2 or by making a least square fit to <rf with the 
appropriate number of harmonics. Note that the definitions of (i) zero order (r = 0) cyclic 
autocorrelation, (ii) periodogram, and (iii) zero order harmonic of the variance, coincide with 
those usually adopted for stationary random processes. Thus, even though a cyclostationary 
time series is not stationary, the ordinary spectral analysis can be used for obtaining the zero 
order spectra. Note, however, that cyclostationary random processes provide more spectral 
information about the time series they are associated with due to the existence of cyclic 
spectra with r > 0. 

As an important and practical application, let us consider a time series y t consisting of 
the sum of a stationary random process, n t , and a cyclostationary one X t (i.e. y t — n t + X t ). 
Let the variance of the stationary time series n t be v 2 and its spectral density be £(/). It 
is easy to see that the resulting process is also cyclostationary. If the two processes are 
uncorrelated, then the zero order harmonic Sq of the variance of the combined processes is 
equal to 

E 2 = u 2 + a 2 , (38) 
and the zero order spectrum, G (f), of y t is 

G (f) = S(f) + g (f) . (39) 

The harmonics of the variance as well as the cyclic spectra of y t with r > coincide instead 
with those of X t . In other words, the harmonics of the variance and the cyclic spectra of the 
process y t with r > contain information only about the cyclostationary process X t , and 
are not "contaminated" by the stationary process n t . 



VI. ANALYTIC STUDY OF THE BACKGROUND SIGNAL 

In the case of the ensemble of N WD- WD binaries, the total signal s(t) is given by the 
following sum 

N 

s(t)=J2x GW (t;A l ) , (40) 



23 



where A represents the set (<f> , l,i[), D, (3, X,J\4 c ,uj s ) of 8 parameters characterizing a GW 
signal. Since N is large, we can expect the parameters of the signals to be randomly 
distributed and regard the signal s(t) itself as a random process. Its mean, m(t), and 
its autocorrelation function, C(t',t) can then be calculated by assuming the probability 
distribution of the vector A, P(A), to be the product of the five probability distributions, 
Pi(0 o ), Pz{l), Psi.i'): Pi(D,(3, A), and P^{M. Ci oj s ) (as we did in our numerical simulation of 
the WD- WD background in Section ITTTj) . By assuming the angles <p and ip to be uniformly 
distributed in the interval [0, 2tt), and cos t to be uniformly distributed in the interval [—1, 1], 
we can then perform the integrals over the angles O , ip, and i analytically and obtain the 
following expressions 

m(t) = N [ X GW (t)P(A)dA 
Jv 

jy r 2ir r 2ir p\ 



/ d(j) I dip I d cos l I / X GW (t)P±P 5 dV 4 dV 5 = , (41) 
Jo Jo J -I Jva Jvr, 



C(t',t) = N [ X GW (t')X GW (t)P(\)dA 
Jv 



= / d<f) [ dip I dcost [ [ ^[A{x,t')A*{x,t)e im) -^ t ' )] ]P i P 5 dV 4 dl^) 

16vr^ Jo Jo J-\ J Vi Jv 5 

Note that the mean value m(t) is equal to as a consequence of averaging the antenna 
response over the polarization angle ip. 

In order to gain an analytic insight about the statistical properties of the autocorrelation 
function C(t', t), in what follows we will adopt the zero-order long- wavelength approximation 
of the LISA response X GW (t) obtained by fixing n — in (|7|) and using the expression for the 
complex amplitude given in equation (|8J) . After some long but straightforward algebra, 
the autocorrelation function, C(t',t), can be written in the following form 

C(f , t) = ™N f [ x'h 2 [u 2 (t')u 2 (t) + v 2 (t')v 2 (t) + u 3 (t')u 3 (t) + v 3 (t')v 3 (t) 
<v A Jv 5 



5 

- u 2 (t')u 3 (t) - v 2 (t')v 3 (t) - u 3 (t')u 2 (t) - v 3 (t')v 2 (t)} cos[0(O - <p(t)]P A P 5 dV A 

where x = u s L, and h Q = iM ^ /z ^J 2//3 (units in which the gravitational constant, G, 
and the speed of light, c, are equal to 1). For frequencies less than 1 mHz the Doppler 
modulation in the phase <f>(t) can be neglected making (p(t) ~ u s t. If we now introduce a 
new time variable r — t' — t and define B(t, r) = C(t + r, t), we have 

B(t, r) = / VM cos(^r) duo s V B r (r) e irnt , (44) 
Jo r= _ 8 
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where 



V(u 8 ) = ^ (oOsL)^ 3 [ (V2M c ) 10/3 Pr o (M c ,uj s ) dM c , (45) 

b J Ma 



and 



D 2 

The functions b r entering into equation (pftjj) are equal to 



B r (r) = I b r (nr) ^ D ^ ' dV 4 . (46) 



V4 



b = V 2 + Ul + (V 2 + U 2 ) cos(fir) + U\ cos(2ftr) + (Vf + U 2 ) cos(3ftr) 

+ (V 4 2 + Ul) cos(4fir) + (V 2 - Ul) cos(4 7o ) , 
61 = e i(Qr/2 - 5o) [(f/ t/i + 7 Fi) cos(ftr/2) + U X U 2 cos(3fir/2) + 

[/ 2 C/ 3 cos(5fir/2) + (C/3C/4 + V^T4) cos(7fir/2) + (-EWi + cos(f2r/2)e* 
b 2 = e i{nT - 2So) [U U 2 cos(ftr) + U 2 U 4 cos(3ftr) + (V^ + U X U Z ) cos(2fir) + 

(-C/ 2 /2 + V 2 /2 - U U 2 cos(ftr))e i4 ^] , 

63 = e i3 ^ T/2 - 5o) [(^o^3 + V r oV 3 )cos(3fir/2) + (f/ 1 f/4 + V r 1 V r 4)cos(5fir/2) + 

((- Wa + V V 3 ) cos(3ftr/2) - C/iZ7 2 cos(fir/2))e i47 °] , 

6 4 = e i2(nr ^ 2,5o) [(?7o^4 + V V 4 ) cos(2fir) + ((V^ - UiU 3 ) cos(fir) + 

(^y - U U 4 ) cos(2fir) - f/ 2 2 /2)e i470 ] , 

b 5 = e i5(nr/2- < 5 )+i4 7 o[(y i y 4 _ cos ( 3 fi r /2) - C/ 2 [/ 3 CO S (Ot/2)] , 

b 6 = e l{3nT - 65o+ ^ o) [-U 2 U 4 cos(ttT) -U 2 /2 + V 2 /2] , 

b 8 = e l4{nT - 2So+ ^[-Ul/2 + V 4 2 /2] , (47) 

where 5o = A — r] , 70 = A — 770 — £0, and the functions Ui, Vi are give in equations (31-39) of 
It is easy to see that the autocorrelation B(t, r) is periodic in t with period one year for 
a fixed r, making it a cyclostationary random process. Note that, if the ecliptic longitude 
A is uniformly distributed, all the coefficients b r given in Eq. (}4Tj) vanish for r > 0, and the 
random process s(t) becomes stationary as the autocorrelation C(t',t) now depends on the 
time difference t' — t. 

The non-stationarity of the WD- WD background was first pointed out by Giampieri and 



Polnarev 



24j under the assumption of sources distributed anisotropically, and they also 



obtained the Fourier expansion of the sample variance and calculated the Fourier coefficient 
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for simplified WD- WD binary distributions in the Galactic disc. What was however not 
realized in their work is that this non-stationary random process is actually cyclostationary, 
i.e. there exists cyclic spectra that can in principle allow us to infer more information about 
the WD- WD background than one could obtain by just estimating the zero-order spectrum. 

If we now set r = in Eq. 021 we obtain the Fourier expansion of the variance cr 2 (t) of 
the cyclostationary process 

8 

a 2 (t) = B(t,0) = B koe lkm , (48) 



where 



B k0 = V I b k0 PA{D ^ X) dV, , (49) 



n 2 



with V = ± J °°V(uj s ) du s , and 



AX) 



u 2 + u 2 + ul + u 2 + ul + V 2 + V 2 + V 2 + V 2 + 



(V 2 -U 2 )cos(A l0 ), (50) 
6 10 = e- i(5o (tWi + UJJ* + U 2 U 3 + U 3 U 4 + V V 1 + V 3 V i + 

(-Uo^ + VoV^ ) , (51) 
b 20 = e- i25o (U U 2 + UiU 3 + U 2 U 4 + V t V 3 + 

{-U 2 /2 + V 2 /2-U U 2 )e^ ) (52) 

b 30 = e- m °(U U 3 + U X U 4 + V V 3 + V X V A + (53) 

{-U Q U 3 -U X U 2 + V V 3 )e^°) 

&4o = e- ii5o (U U 4 + V V 4 + (54) 

(-U U 4 - U X U 3 + V V 4 + VtV 3 - U 2 /2)e^°) , 

ho = e'^-^i-UiUi - U 2 U 3 + V 1 V 4 ) , (55) 

b m = e^°- 6<5 °)(-[/ 2 £/ 4 - U 2 /2 + V 2 /2) , (56) 

b 70 = e^-^i-UM + V 3 V 4 ) , (57) 

bso = e^- 8So \-U 2 j2 + Vi/2) . (58) 

If we assume the function V(uj s ) to change very little over a frequency bin, or equivalently 
choose t to be such that Qt <C 1 , we can then approximate the functions b r with the functions 
b r o- Under this approximation the cyclic spectra of the process s(t) can be shown to reduce 
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to the following expression 

9r(u s ) = -V(uj s )B Qr . (59) 

Thus under the above approximations the cyclic spectra are determined by one function 
of the Fourier frequency, and by the coefficients of the Fourier decomposition of the cyclic 
variance. Note that this simplified representation of the cyclic spectra will not be valid 
if there are additional correlations between the parameters of the binary population. For 
example, if the chirp masses or the frequencies of the radiation emitted by the binaries are 
correlated with the positions of the binaries themselves in the Galactic disc, then the cyclic 
spectra will display a different frequency dependence from that implied by equation (|59)l . In 
general we can expect the direct measurements of the cyclic spectra from the LISA data to 
allow us to infer properties of the distribution of the parameters characterizing the WD- WD 
population. In other words, by analyzing the 17 real and independent cyclic spectra we 
should be able to derive more information about the WD- WD binary population than we 
would have by simply looking at the ordinary spectrum. 



VII. DATA ANALYSIS OF THE BACKGROUND SIGNAL 

We have numerically implemented the methods outlined in Section IS and applied them 
to our simulated WD- WD background signal. A comparison of the results of our simulation 
of the detached WD- WD background with the calculation of the background by Hils and 



Bender [19. |29| is shown in Figure (jHJ). We find that the amplitude of the background from 
our simulation is a factor of more than 2 smaller than that of Hils and Bender. The level of 
the WD- WD background is determined by the number of such systems in the Galaxy. We 
estimate that our number WD- WD binaries should be correct within a factor 5 and thus the 
amplitude of the background should be right within a factor of \/5. In Figure © we have 
plotted the two backgrounds against the LISA spectral density and we have also included the 
LISA sensitivity curve. The latter is obtained by dividing the instrumental noise spectral 
density by the detector GW transfer function averaged over isotropically distributed and 
randomly polarized signals. In the zero-order long wavelength approximation this averaged 



transfer function is equal to a/3/20. 

Our analysis was applied to 3 years of LISA X data consisting of a coherent superposition 
of signals emitted by detached WD- WD binaries, by semi-detached binaries (AM CVn sys- 



27 



Spectrum of the LISA noise vs. WDWD background noise 
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FIG. 9: Comparison of detached WD- WD background obtained from binary population synthesis 
simulation ( P, ^|) with the WD- WD background calculated by Hils and Bender The 
amplitude spectral density of the LISA instrumental noise and the LISA sensitivity curve are 
drawn for comparison. All spectral densities are one-sided. 

terns), and of simulated instrumental noise. The noise was numerically generated by using 
the spectral density of the TDI X observable given in Q|. In addition a 1 mHz low-pass 
filter was applied to our data set in order to focus our analyses to the frequency region in 
which the WD- WD stochastic background is expected to be dominant. 

The results of the Fourier analysis of the sample variance of the background signal are 
shown in Figures (|1U|) and (J 11)) . The top panel of Figure (|1U|) shows the sample variance of the 
simulated data for which the variances were estimated over a period of 1 week; periodicity is 
clearly visible. The bottom panel instead shows the Fourier analysis of the sample variance 
for which we have removed the mean from the sample variance time series. The vertical 
lines correspond to multiples of 1 year; two harmonics can clearly be distinguished from 
noise. The other peaks of the spectrum that fall roughly half way between the multiples 
of 1/year frequency, are from the rectangular window inherent to the finite time series. In 
Figure (fTTjl we present the least square fit of 8 harmonics to our 3 years of simulated X 
data. The number 8 comes from our theoretical predictions of the number of harmonics 
obtained in Section IVTl fsee Eq. (}4*4*j) b We have calculated the magnitude of the harmonics 
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WD-WD background - 3 years of data 




Fourier analysis of the variance 




FIG. 10: Top panel: The sample variance of the simulated WD-WD background observed by LISA. 
The data includes two populations of WD-WD binaries, detached and semi-detached, which are 
added to the LISA instrumental noise. The data is passed through a low-pass filter with a cut-off 
frequency of 1 mHz. Bottom panel: Fourier analysis of the sample variance. Two harmonics are 
clearly resolved. 



and obtained the residuals. The results from the least square fit agree very well with those 
obtained via Fourier analysis (see also Figure (|12p). The magnitudes of the first and second 
harmonics resolved by Fourier analysis, for instance, agree with the corresponding least 
square fit estimates within a few hundredth of a percent. 

It is useful to compare the results of our numerical analysis against the analytic calcu- 
lations of Giampieri and Polnarev 24]. Their analytic expressions for the harmonics of the 
variance of a background due to binary systems distributed in the galactic disc are given in 
Eq. (42) and shown in Figure (4) of 24]. Our estimation roughly matches theirs in that 
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FIG. 11: Top panel: The sample variance of the WD- WD background data and the least square fit 
to it using 8 harmonics (small circles). Middle panel: magnitude of the harmonics obtained from 
the least square fit. Bottom panel: residual error between the fit and the data. 



the 0th order harmonics is dominant and the first two harmonics have more power than the 
remaining ones. Our estimate of the power in the second harmonic, however, is larger than 
that in the first one, whereas they find the opposite. We attribute this difference to their 
use of a Gaussian distribution of sources in the Galactic disc rather than the exponential 
that we adopted from Q]. Comparison between these two results suggests that it should be 
possible to infer the distribution of WD-WD binaries in our Galaxy by properly analyzing 
the harmonics of the variance of the galactic background measured by LISA. How this can 
be done will be the subject of a future work. 
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FIG. 12: Comparison between the estimated power in the harmonics obtained via (i) Fourier 
analysis, (ii) least square fit, and (iii) numerical calculation based on Eq. 1491 The blue line is the 
power spectrum of the variance of the data, red vertical lines are obtained from least square fit 
and the black vertical lines are from the numerical calculation. 



In order to validate our simulation and data analysis method we have compared the 
results of our estimation of the power in the harmonics of the variance against the explicit 
analytic calculation. To estimate the powers we have used Eq. and we have evaluated 
the integrals by numerical and Monte Carlo methods. In the numerical calculation of the 
harmonics we have limited our analysis to the population of detached WD- WD binaries. 
Thus in order to make the comparison meaningful we have performed Fourier analysis and 
least square fit of the time series consisting only of simulated detached WD-WD binaries 
(without semi-detached ones and LISA instrumental noise). The results of the comparison 
are given in Figure (J12)) . We see that for the Oth order harmonic and the first two harmonics 
the agreement is very good. For higher order harmonics there are large discrepancies between 
the numerical calculation and estimation by the least square fit, while by using the Fourier 
transform method, we cannot even resolve higher harmonics in our 3-year data set. We 
conclude that only the two first harmonics can be extracted reliably from the data. We 
also observe a very good agreement between the Fourier and the least square method. As 
a next step in our analysis, we have estimated the cyclic spectra of the simulated WD-WD 
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Cyclic spectra of WD-WD background 
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FIG. 13: The main (k = 0) spectrum of the simulated WD-WD background signal (red), and the 
8 cyclic spectra (magenta) estimated from the simulated data are shown. The spectral density of 
the LISA instrumental noise (black) is shown for reference. 

background signal. In Figure (fT3*j) we have shown the cyclic spectra estimated from the 
data. We have also plotted the spectrum of the LISA instrumental noise and the main 
spectrum (k = 0) estimated from the simulation. We find that the main spectrum and two 
cyclic spectra for k = 1 and k = 2 have the largest magnitude and, over some frequency 
range, they lie above the LISA instrumental noise. The remaining spectra are an order of 
magnitude smaller and are very noisy. We also see that all the cyclic spectra have roughly 
the same slope. This is predicted by our analytic calculations in Section EU and it follows 
from the assumed independence between the location of the binaries in the Galaxy (D, A, (3) 
and their frequencies and chirp masses (u s ,Ai c )- We also find the magnitude of the 2nd 
cyclic spectrum to be higher than the first, similarly to what we had for the harmonics of 
the variance. Note that we estimated the spectra from the time series consisting of the 
WD-WD background added to the LISA instrumental noise. Like the analysis we did for 
the variance, we have also compared the estimates of the cyclic spectra from our simulation 
against those obtained via numerical calculation of the equations derived in Section IVTl The 
corresponding results are presented in Figures (fPfj) and (fT5j) . where it is shown that the 
agreement between the two is quite good. 
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FIG. 14: Estimated main {k = 0) spectrum of the WD- WD background (red) against the calculated 
spectrum (black). The LISA spectral density curve (blue) is shown for comparison. The 0th order 
spectrum contains the LISA instrumental and hence it differs from the spectrum given in Figure 

Our analysis has shown that the LISA data will allow us to compute 17 independent 
cyclic spectra (the 8 complex cyclic spectra g r (f), t = 1, 2, ...8 and the real spectrum go(f)) 
of the WD-WD galactic background, 5 of which can be expected to be measured reliably. 
We have also shown that by performing generalized spectral analysis of the LISA data we 
will be able to derive more information about the WD-WD binary population (properties 
of the distribution of its parameters) than we would have by only looking at the ordinary 
9o(f) spectrum. 
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FIG. 15: Estimated cyclic spectra (black) against the calculated spectra (red). The LISA spectral 
density curve (blue) is shown for comparison. 
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